close , clear, clc;
%给定参数
Nt = 80; 
p1 = [100, 100, 100, 100, 100, 100, 100, 100]; 
p2 = [50, 50, 50, 50, 50, 50, 50, 50]; 
p3 = [100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000];
h1 = 100; 
h2 = [h1/5, h1/3, h1/2, h1, h1*2, h1*3, h1*5, h1*24];
u2 = p2 ./ p1; u3 = p3 ./ p2; 
v2 = h2 ./ h1;
i = sqrt(-1);
%进行计算
for ix = 1:8
    for  iy = 1: Nt
        B(ix,iy) = 10 ^ ((iy+10) / 10) / h1;
        C(ix,iy) = ( coth( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) + acoth( sqrt(u2(ix)) * ...      %双曲余切
            coth( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) * v2(ix) / sqrt(u2(ix)) + acoth( sqrt( u3(ix)/u2(ix)))))))^2;
        D(ix,iy) = ( tanh( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) + atanh( sqrt(u2(ix)) * ...      %双曲正切
            tanh( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) * v2(ix) / sqrt(u2(ix)) + atanh( sqrt( u3(ix)/u2(ix)))))))^2;
    end
    Re1(ix,:) = real(C(ix,:));          %实部
    Im1(ix,:) = imag(D(ix,:));          %虚部
end
%%
%进行画图
figure(1)

loglog(B(1,:),Re1(1,:))
title("H型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")
hold on

loglog(B(2,:),Re1(2,:))
title("H型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

hold on
loglog(B(3,:),Re1(3,:))
title("K型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

hold on
loglog(B(4,:),Re1(4,:))
title("K型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

hold on
loglog(B(5,:),Re1(5,:))
title("A型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

hold on
loglog(B(6,:),Re1(6,:))
title("A型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

hold on
loglog(B(7,:),Re1(7,:))
title("Q型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

hold on
loglog(B(8,:),Re1(8,:))
title("Q型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

hold off
